MODULE VSURFSAD_MOD
CONTAINS
SUBROUTINE VSURFSAD(KIDIA,KFDIA,KLON,KLEVS,KTILE,&
 & KTVL,KTVH,&
 & PLAIL, PLAIH,&
 & PTMLEV5 ,PQMLEV5  ,PAPHMS5 ,&
 & PTSKM1M5,PWSAM1M5 ,PTSAM1M5,KSOTY,&
 & PSRFD5  ,PRAQ5    ,&
 & PSSDP2, PSSDP3, YDCST   ,YDVEG    ,YDSOIL  ,&
 & PQSAM5  ,PQS5    ,PDQS5    ,&
 & PWETB5  ,PCPTS5   ,PWETL5  , PWETH5, PWETHS5, &
 & PTMLEV  ,PQMLEV   ,PAPHMS  ,&
 & PTSKM1M ,PWSAM1M  ,PTSAM1M ,&
 & PSRFD   ,PRAQ     ,PQSAM   ,&
 & PQS     ,PDQS     ,&
 & PWETB   ,PCPTS    ,PWETL   , PWETH , PWETHS)

USE PARKIND1 , ONLY : JPIM, JPRB
USE YOMHOOK  , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_THF  , ONLY : R4LES, R5LES, R2ES, R4IES, R3LES, R3IES, R5IES, RVTMP2
USE YOS_CST  , ONLY : TCST
USE YOS_VEG  , ONLY : TVEG
USE YOS_SOIL , ONLY : TSOIL
USE YOMSURF_SSDP_MOD

! (C) Copyright 2006- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!     ------------------------------------------------------------------

!**   *VSURFSAD* - PREPARES SURFACE BOUNDARY CONDITION FOR T AND Q
!                   (Adjoint)

!  AUTHOR:
!    M. Janiskova       ECMWF   March 2006
!
!  Modified:
!    G. Balsamo         ECMWF   January 2007 Add soil type
!    M. Janiskova       ECMWF   August 2007  Cleanning for soil type
!    S. Boussetta/G.Balsamo May 2009    Add lai
!    M. Janiskova           August 2011 Security for lai
!    I. Ayan-Miguez July 2023 Added PSSDP2 and PSSDP3 objects for spatially distributed parameters
!     PURPOSE
!     -------

!     PREPARE SURFACE BOUNDARY CONDITION FOR Q AND T, E.G. FRACTIONAL
!     SURFACE COVER (SNOW AND VEGETATION), SATURATION SPECIFIC HUMIDITY
!     AT THE SURFACE, RELATIVE HUMIDITY OVER BARE LAND AND THE STOMATAL
!     RESISTANCE.

!     INTERFACE
!     ---------

!     *VSURFSAD* IS CALLED BY *SURFEXCDRIVERSAd*

!     INPUT PARAMETERS (INTEGER):

!     *KIDIA*        START POINT
!     *KFDIA*        END POINT
!     *KLON*         NUMBER OF GRID POINTS PER PACKET
!     *KLEVS*        NUMBER OF SOIL LAYERS
!     *KTILE*        TILE INDEX
!     *KTVL*         VEGETATION TYPE FOR LOW VEGETATION FRACTION
!     *KTVH*         VEGETATION TYPE FOR HIGH VEGETATION FRACTION
!     *PLAIL*        LAI OF LOW VEGETATION
!     *PLAIH*        LAI OF HIGH VEGETATION
!     *KSOTY*        SOIL TYPE 

!     INPUT PARAMETERS (REAL):
!  Trajectory  Perturbation  Description                               Unit
!  PTMLEV5     PTMLEV        TEMPERATURE AT T-1, lowest model level    K
!  PQMLEV5     PQMLEV        SPECIFIC HUMIDITY AT T-1, lowest model 
!                            level                                     kg/kg
!  PAPHMS5     PAPHMS        PRESSURE AT T-1, surface                  Pa
!  PTSKM1M5    PTSKM1M       SURFACE TEMPERATURE                       K
!  PWSAM1M5    PWSAM1M       SOIL MOISTURE ALL LAYERS                 m**3/m**3
!  PTSAM1M5    PTSAM1M       SOIL TEMPERATURE ALL LAYERS               K
!  PSRFD5      PSRFD         DOWNWARD SHORT WAVE RADIATION FLUX AT 
!                            SURFACE                                   W/m2
!  PRAQ5       PRAQ          PRELIMINARY AERODYNAMIC RESISTANCE

!     OUTPUT PARAMETERS (REAL):
!  Trajectory  Perturbation  Description
!  PQSAM5      PQSAM         SPECIFIC HUMIDITY AT THE SURFACE          kg/kg
!  PQS5        PQS           SATURATION Q AT SURFACE                   kg/kg
!  PDQS5       PDQS          DERIVATIVE OF SATURATION Q-CURVE AT
!                            SURFACE T
!  PWETB5      PWETB         BARE SOIL RESISTANCE
!  PCPTS5      PCPTS         DRY STATIC ENRGY AT SURFACE               J/kg
!  PWETL5      PWETL         CANOPY RESISTANCE LOW VEGETATION
!  PWETH5      PWETH         CANOPY RESISTANCE HIGH VEGETATION, 
!                            SNOW FREE
!  PWETHS5     PWETHS        CANOPY RESISTANCE HIGH VEGETATION 
!                            WITH SNOW

!     METHOD
!     ------

!     SEE DOCUMENTATION

!     ------------------------------------------------------------------

IMPLICIT NONE

INTEGER(KIND=JPIM),INTENT(IN)    :: KLON 
INTEGER(KIND=JPIM),INTENT(IN)    :: KLEVS 
INTEGER(KIND=JPIM),INTENT(IN)    :: KIDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KFDIA 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTILE 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTVL(:) 
INTEGER(KIND=JPIM),INTENT(IN)    :: KTVH(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PLAIL(:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PLAIH(:)
INTEGER(KIND=JPIM),INTENT(IN)    :: KSOTY(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PQMLEV5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PAPHMS5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSKM1M5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PWSAM1M5(:,:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PTSAM1M5(:,:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSRFD5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PRAQ5(:) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSSDP2(:,:)
REAL(KIND=JPRB)   ,INTENT(IN)    :: PSSDP3(:,:,:)
TYPE(TCST)        ,INTENT(IN)    :: YDCST
TYPE(TVEG)        ,INTENT(IN)    :: YDVEG
TYPE(TSOIL)       ,INTENT(IN)    :: YDSOIL
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PQSAM5(:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PQS5(:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PDQS5(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETB5(:) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: PCPTS5(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETL5(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETH5(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETHS5(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PTMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PQMLEV(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PAPHMS(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PTSKM1M(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWSAM1M(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PTSAM1M(:,:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PSRFD(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PRAQ(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PQSAM(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PQS(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PDQS(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETB(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PCPTS(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETL(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETH(:) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: PWETHS(:) 

!*    LOCAL STORAGE
!     ----- -------

REAL(KIND=JPRB) :: ZLIQ(KLON,KLEVS)
REAL(KIND=JPRB) :: ZLIQ5(KLON,KLEVS)
REAL(KIND=JPRB) :: ZF5(KLON,KLEVS),ZLIQP15(KLON,KLEVS),ZLIQP25(KLON,KLEVS)

REAL(KIND=JPRB) :: ZRSMINH(KLON), ZRSMINL(KLON)
REAL(KIND=JPRB) :: ZLAIH(KLON), ZLAIL(KLON)
REAL(KIND=JPRB) :: ZROOT1L(KLON),ZROOT2L(KLON),ZROOT3L(KLON),ZROOT4L(KLON)
REAL(KIND=JPRB) :: ZROOT1H(KLON),ZROOT2H(KLON),ZROOT3H(KLON),ZROOT4H(KLON)
REAL(KIND=JPRB) :: ZHSTRH(KLON), ZHSTRL(KLON)

REAL(KIND=JPRB) :: ZCOR5(KLON), ZWROOTL5(KLON), ZWROOTH5(KLON)
REAL(KIND=JPRB) :: ZF21LP15(KLON), ZF21LP25(KLON), ZF21L5(KLON), ZF2L5(KLON)
REAL(KIND=JPRB) :: ZF21HP15(KLON), ZF21HP25(KLON), ZF2H5(KLON), ZF21H5(KLON)
REAL(KIND=JPRB) :: ZF21BP15(KLON), ZF21BP25(KLON), ZF21B5(KLON), ZF2B5(KLON)
REAL(KIND=JPRB) :: ZSRFL5(KLON), ZQSAIR5(KLON), ZQSAIR55(KLON)
REAL(KIND=JPRB) :: ZF1LPI5(KLON), ZF1LP5(KLON)
REAL(KIND=JPRB) :: ZF1L5(KLON), ZF3H5(KLON), ZF3L5(KLON), ZF1H5(KLON)
REAL(KIND=JPRB) :: ZWETL5(KLON), ZWETLI5(KLON)
REAL(KIND=JPRB) :: ZWETH5(KLON), ZWETHI5(KLON), ZF2BI5(KLON)
REAL(KIND=JPRB) :: Z3S5(KLON), Z4S5(KLON), Z5S5(KLON), Z6S5(KLON)
REAL(KIND=JPRB) :: ZDIV15(KLON), ZDIV25(KLON), ZDIV35(KLON)
REAL(KIND=JPRB) :: ZDIV45(KLON), ZDIV55(KLON), ZDIV65(KLON)
REAL(KIND=JPRB) :: ZDIV75(KLON), ZDIV85(KLON), ZDIV95(KLON)
REAL(KIND=JPRB) :: ZCOR45(KLON), ZF3L55(KLON), ZF3H55(KLON)
REAL(KIND=JPRB) :: ZPQS5(KLON)

INTEGER(KIND=JPIM) :: JK, JL, JS

REAL(KIND=JPRB) ::   &
 & ZCOR, ZEPSF3, ZF, ZF1H, ZF1L, ZF2H, ZF2L, ZF2B, ZF21H, ZF21L, ZF21B, ZF3H, ZF3L,   &
 & ZQSAIR, ZRVA, ZRVB, ZSRFL, ZWROOTH, ZWROOTL
REAL(KIND=JPRB) :: ZQWEVAP, ZWPWP

REAL(KIND=JPRB) :: Z3S, Z4S, ZLIQP1, ZLIQP2, ZF2LP1, ZF2LP2,ZF21LP1, ZF21LP2
REAL(KIND=JPRB) :: ZF2HP1, ZF2HP2, ZF21HP1, ZF21HP2, ZF2BP1, ZF2BP2, ZF21BP1, ZF21BP2, ZF1LP
REAL(KIND=JPRB) :: Z5S, Z6S, ZWETL, ZWETH

REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!     ------------------------------------------------------------------

!*       1.     INITIALIZE CONSTANTS
!               ---------- ----------

IF (LHOOK) CALL DR_HOOK('VSURFSAD_MOD:VSURFSAD',0,ZHOOK_HANDLE)
ASSOCIATE(RCPD=>YDCST%RCPD, RETV=>YDCST%RETV, RLSTT=>YDCST%RLSTT, &
 & RLVTT=>YDCST%RLVTT, RTT=>YDCST%RTT, &
 & LEVGEN=>YDSOIL%LEVGEN, RQWEVAP=>YDSOIL%RQWEVAP, &
 & RTF1=>YDSOIL%RTF1, RTF2=>YDSOIL%RTF2, RTF3=>YDSOIL%RTF3, RTF4=>YDSOIL%RTF4, &
 & RWCAP=>YDSOIL%RWCAP, RWCAPM3D=>PSSDP3(:,:,SSDP3D_ID%NRWCAPM3D), RWPWP=>YDSOIL%RWPWP, &
 & RWPWPM3D=>PSSDP3(:,:,SSDP3D_ID%NRWPWPM3D), RCEPSW=>YDVEG%RCEPSW, RVLAI=>YDVEG%RVLAI, &
 & RVROOTSAL3D=>PSSDP3(:,:,SSDP3D_ID%NRVROOTSAL3D), RVROOTSAH3D=>PSSDP3(:,:,SSDP3D_ID%NRVROOTSAH3D), &
 & RVRSMINL2D=>PSSDP2(:,SSDP2D_ID%NRVRSMINL2D),  RVRSMINH2D=>PSSDP2(:,SSDP2D_ID%NRVRSMINH2D), &
 & RWRESTM3D=>PSSDP3(:,:,SSDP3D_ID%NRWRESTM3D), RVRSMINB2D=>PSSDP2(:,SSDP2D_ID%NRVRSMINB2D), &
 & RRSF1A=>YDSOIL%RRSF1A, RRSF1B=>YDSOIL%RRSF1B, RRSF1C=>YDSOIL%RRSF1C, &
 & RVHSTRL2D=>PSSDP2(:,SSDP2D_ID%NRVHSTRL2D), RVHSTRH2D=>PSSDP2(:,SSDP2D_ID%NRVHSTRH2D))

ZRVA=5000._JPRB
ZRVB=10._JPRB
ZEPSF3=0.00001_JPRB ! security value for exponential sat-deficit dependence

!     ------------------------------------------------------------------

!          2.    PREPARE SURFACE BOUNDARY CONDITION
!                ------- ------- -------- ---------

!*         2.1   RELATIVE HUMIDITY OVER THE BARE LAND PART

!                BARE SOIL RESISTANCE IS COMPUTED FOR KTILE=4

!*         2.2   SATURATION PARAMETERS,

DO JL=KIDIA,KFDIA
  IF (PTSKM1M5(JL) > RTT) THEN
    ZDIV15(JL) = 1.0_JPRB/(PTSKM1M5(JL)-R4LES)
    Z3S5(JL) = R3LES*(PTSKM1M5(JL)-RTT)*ZDIV15(JL)
  ELSE
    ZDIV15(JL) = 1.0_JPRB/(PTSKM1M5(JL)-R4IES)
    Z3S5(JL) = R3IES*(PTSKM1M5(JL)-RTT)*ZDIV15(JL)
  ENDIF
  Z4S5(JL) = EXP(Z3S5(JL))
  ZDIV25(JL) = 1.0_JPRB/PAPHMS5(JL)
  ZPQS5(JL) = (R2ES*Z4S5(JL))*ZDIV25(JL)
  ZCOR5(JL) = 1.0_JPRB/(1.0_JPRB-RETV*ZPQS5(JL))
  PQS5(JL) = ZPQS5(JL)*ZCOR5(JL)

  IF (PTSKM1M5(JL) > RTT) THEN
    ZDIV35(JL) = 1.0_JPRB/(PTSKM1M5(JL)-R4LES)
    PDQS5(JL) = R5LES*PQS5(JL)*ZCOR5(JL)*ZDIV35(JL)*ZDIV35(JL)
  ELSE
    ZDIV35(JL) = 1.0_JPRB/(PTSKM1M5(JL)-R4IES)
    PDQS5(JL) = R5IES*PQS5(JL)*ZCOR5(JL)*ZDIV35(JL)*ZDIV35(JL)
  ENDIF
ENDDO

!*         2.3   DEFINITION OF THE STOMATAL RESISTANCE AND BARE SOIL RES
!*               DOES WORK FOR TYPE 4, 6 AND 8 WHEN ROUTINE IS CALLED FOR 
!*               TYPE 4

IF (KTILE  ==  4) THEN

!                Compute first liquid fraction of soil water to 
!                be used later in stress functions
!          CONTRIBUTION TO APPARENT ENERGY, TAKING INTO ACCOUNT
!          FREEZING/MELTING OF SOIL WATER.

  DO JK=1,KLEVS
    DO JL=KIDIA,KFDIA
      IF (PTSAM1M5(JL,JK) < RTF1.AND.PTSAM1M5(JL,JK) > RTF2) THEN
        ZF5(JL,JK) = 0.5_JPRB*(1.0_JPRB-SIN(RTF4*(PTSAM1M5(JL,JK)-RTF3)))
      ELSEIF (PTSAM1M5(JL,JK) <= RTF2) THEN
        ZF5(JL,JK) = 1.0_JPRB
      ELSE
        ZF5(JL,JK) = 0.0_JPRB
      ENDIF

      ZLIQP15(JL,JK) = PWSAM1M5(JL,JK)*(1.0_JPRB-ZF5(JL,JK))
      IF (LEVGEN) THEN
        JS=KSOTY(JL)
        IF (ZLIQP15(JL,JK) < RWCAPM3D(JL,JK)) THEN
          ZLIQP25(JL,JK) = ZLIQP15(JL,JK)
        ELSE
          ZLIQP25(JL,JK) = RWCAPM3D(JL,JK)
        ENDIF
        IF (ZLIQP25(JL,JK) > RWPWPM3D(JL,JK)) THEN
          ZLIQ5(JL,JK) = ZLIQP25(JL,JK)
        ELSE
          ZLIQ5(JL,JK) = RWPWPM3D(JL,JK)
        ENDIF
      ELSE
        IF (ZLIQP15(JL,JK) < RWCAP) THEN
          ZLIQP25(JL,JK) = ZLIQP15(JL,JK)
        ELSE
          ZLIQP25(JL,JK) = RWCAP
        ENDIF
        IF (ZLIQP25(JL,JK) > RWPWP) THEN
          ZLIQ5(JL,JK) = ZLIQP25(JL,JK)
        ELSE
          ZLIQ5(JL,JK) = RWPWP
        ENDIF
      ENDIF
    ENDDO
  ENDDO

  DO JL=KIDIA,KFDIA
!           minimal stomatal resistance : ZRSMIN
    ZRSMINL(JL) = RVRSMINL2D(JL)
    ZRSMINH(JL) = RVRSMINH2D(JL)

!           leaf area index  : ZLAI
    ZLAIL(JL) = PLAIL(JL)
    ZLAIH(JL) = PLAIH(JL)

!           soil moisture stress function : F2
    ZROOT1L(JL) = RVROOTSAL3D(JL,1)
    ZROOT2L(JL) = RVROOTSAL3D(JL,2)
    ZROOT3L(JL) = RVROOTSAL3D(JL,3)
    ZROOT4L(JL) = RVROOTSAL3D(JL,4)
    ZROOT1H(JL) = RVROOTSAH3D(JL,1)
    ZROOT2H(JL) = RVROOTSAH3D(JL,2)
    ZROOT3H(JL) = RVROOTSAH3D(JL,3)
    ZROOT4H(JL) = RVROOTSAH3D(JL,4)

    ZWROOTL5(JL) = ZLIQ5(JL,1)*ZROOT1L(JL)+ZLIQ5(JL,2)*ZROOT2L(JL) &
     &       + ZLIQ5(JL,3)*ZROOT3L(JL)+ZLIQ5(JL,4)*ZROOT4L(JL)
    ZWROOTH5(JL) = ZLIQ5(JL,1)*ZROOT1H(JL)+ZLIQ5(JL,2)*ZROOT2H(JL) &
     &       + ZLIQ5(JL,3)*ZROOT3H(JL)+ZLIQ5(JL,4)*ZROOT4H(JL)

    IF (LEVGEN) THEN
      JS=KSOTY(JL)
      ZWPWP=RWPWPM3D(JL,1)
      IF (JS >= 1) THEN
        ZQWEVAP=1._JPRB/(RWCAPM3D(JL,1_JPIM)-ZWPWP)
      ELSE
        ZQWEVAP = 0.0_JPRB
      ENDIF
    ELSE
      ZWPWP=RWPWP
      ZQWEVAP=RQWEVAP
    ENDIF

    IF (PTSAM1M5(JL,1) <= RTT ) THEN
      ZF2L5(JL)=RCEPSW
      ZF2H5(JL)=RCEPSW
    ELSE
      ZF21LP15(JL) = (ZWROOTL5(JL)-ZWPWP)*ZQWEVAP
      IF (ZF21LP15(JL) < 1.0_JPRB) THEN
        ZF21LP25(JL) = ZF21LP15(JL)
      ELSE
        ZF21LP25(JL) = 1.0_JPRB
      ENDIF
      IF (ZF21LP25(JL) > RCEPSW) THEN
        ZF21L5(JL) = ZF21LP25(JL)
      ELSE
        ZF21L5(JL) = RCEPSW
      ENDIF

      ZF2L5(JL) = 2._JPRB*ZF21L5(JL)-(ZF21L5(JL)*ZF21L5(JL)) 


      ZF21HP15(JL) = (ZWROOTH5(JL)-ZWPWP)*ZQWEVAP
      IF (ZF21HP15(JL) < 1.0_JPRB) THEN
        ZF21HP25(JL) = ZF21HP15(JL)
      ELSE
        ZF21HP25(JL) = 1.0_JPRB
      ENDIF
      IF (ZF21HP25(JL) > RCEPSW) THEN
        ZF21H5(JL) = ZF21HP25(JL)
      ELSE
        ZF21H5(JL) = RCEPSW
      ENDIF

      ZF2H5(JL) = 2._JPRB*ZF21H5(JL)-(ZF21H5(JL)*ZF21H5(JL)) 

    ENDIF

    ZF21BP15(JL) = (ZLIQ5(JL,1)-ZWPWP)*ZQWEVAP
    IF (ZF21BP15(JL) < 1.0_JPRB) THEN
      ZF21BP25(JL) = ZF21BP15(JL)
    ELSE
      ZF21BP25(JL) = 1.0_JPRB
    ENDIF
    IF (ZF21BP25(JL) > RCEPSW) THEN
      ZF21B5(JL) = ZF21BP25(JL)
    ELSE
      ZF21B5(JL) = RCEPSW
    ENDIF

    ZF2B5(JL) = 2._JPRB*ZF21B5(JL)-(ZF21B5(JL)*ZF21B5(JL)) 

!           radiation stress function (proposed by Alan Betts): ZF1 

    ZSRFL5(JL) = PSRFD5(JL)/RRSF1B
    ZF1LPI5(JL) = 1.0_JPRB/(ZSRFL5(JL)+RRSF1C)
    ZF1LP5(JL) = RRSF1A*(1.+ZSRFL5(JL))*ZF1LPI5(JL)
    IF (ZF1LP5(JL) > 1.0_JPRB) THEN
      ZF1L5(JL) = 1.0_JPRB/ZF1LP5(JL)
    ELSE
      ZF1L5(JL) = 1.0_JPRB
    ENDIF
    ZF1H5(JL) = ZF1L5(JL)

!           atmospheric moisture deficit stress function : F3
    ZHSTRL(JL)=RVHSTRL2D(JL)
    ZHSTRH(JL)=RVHSTRH2D(JL)

    IF (PTMLEV5(JL) > RTT) THEN
      ZDIV45(JL) = 1.0_JPRB/(PTMLEV5(JL)-R4LES)
      Z5S5(JL) = R3LES*(PTMLEV5(JL)-RTT)*ZDIV45(JL)
    ELSE
      ZDIV45(JL) = 1.0_JPRB/(PTMLEV5(JL)-R4IES)
      Z5S5(JL) = R3IES*(PTMLEV5(JL)-RTT)*ZDIV45(JL)
    ENDIF
    Z6S5(JL) = EXP(Z5S5(JL))
    ZDIV55(JL) = 1.0_JPRB/PAPHMS5(JL)
    ZQSAIR55(JL) = (R2ES*Z6S5(JL))*ZDIV55(JL)
    ZCOR45(JL) = 1.0_JPRB/(1.0_JPRB-RETV*ZQSAIR55(JL))
    ZQSAIR5(JL) = ZQSAIR55(JL)*ZCOR45(JL)

    ZF3L5(JL) = EXP(-ZHSTRL(JL)*(ZQSAIR5(JL)-PQMLEV5(JL)))
    ZF3H5(JL) = EXP(-ZHSTRH(JL)*(ZQSAIR5(JL)-PQMLEV5(JL)))
    ZF3L55(JL) = ZF3L5(JL)
    ZF3H55(JL) = ZF3H5(JL)
    IF (ZF3L55(JL) > 1.0_JPRB) THEN
      ZF3L5(JL) = 1.0_JPRB
    ENDIF
    IF (ZF3L55(JL) < ZEPSF3) THEN
      ZF3L5(JL) = ZEPSF3
    ENDIF
    IF (ZF3H55(JL) > 1.0_JPRB) THEN
      ZF3H5(JL) = 1.0_JPRB
    ENDIF
    IF (ZF3H55(JL) < ZEPSF3) THEN
      ZF3H5(JL) = ZEPSF3
    ENDIF

    IF (ZLAIL(JL) /= 0.0_JPRB) THEN
      ZWETL5(JL) = ZLAIL(JL)*ZF1L5(JL)*ZF2L5(JL)*ZF3L5(JL)
      ZWETLI5(JL) = 1.0_JPRB/ZWETL5(JL)
      PWETL5(JL) = ZRSMINL(JL)*ZWETLI5(JL)
    ELSE
      PWETL5(JL) = 1.0E+6_JPRB
    ENDIF
    IF (ZLAIH(JL) /= 0.0_JPRB) THEN
      ZWETH5(JL) = ZLAIH(JL)*ZF1H5(JL)*ZF2H5(JL)*ZF3H5(JL)
      ZWETHI5(JL) = 1.0_JPRB/ZWETH5(JL)
      PWETH5(JL) = ZRSMINH(JL)*ZWETHI5(JL)
    ELSE
      PWETH5(JL) = 1.0E+6_JPRB  
    ENDIF

    PWETHS5(JL) = PWETH5(JL)
    ZF2BI5(JL) = 1.0_JPRB/ZF2B5(JL)
    PWETB5(JL) = RVRSMINB2D(JL)*ZF2BI5(JL)
  ENDDO
ENDIF

IF (KTILE == 4) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV5(JL) > PQS5(JL)) THEN
      PWETL5(JL) = 0.0_JPRB
    ENDIF
  ENDDO
ELSEIF (KTILE == 6) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV5(JL) > PQS5(JL)) THEN
      PWETH5(JL) = 0.0_JPRB
    ENDIF
  ENDDO
ELSEIF (KTILE == 7) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV5(JL) > PQS5(JL)) THEN
      PWETHS5(JL) = 0.0_JPRB
    ENDIF
  ENDDO
ELSEIF (KTILE == 8) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV5(JL) > PQS5(JL)) THEN
      PWETB5(JL) = 0.0_JPRB
    ENDIF
  ENDDO
ENDIF

!*         2.4   APPARENT SURFACE HUMIDITY

IF (KTILE  ==  1.OR. KTILE  ==  2.OR. KTILE  ==  3.OR. KTILE  ==  5) THEN
  DO JL=KIDIA,KFDIA
    PQSAM5(JL) = PQS5(JL)
  ENDDO
ELSEIF (KTILE  ==  8) THEN
  DO JL=KIDIA,KFDIA
    ZDIV65(JL) = 1.0_JPRB/(PWETB5(JL)+PRAQ5(JL))
    PQSAM5(JL) = PQS5(JL)+(PQMLEV5(JL)-PQS5(JL))*PWETB5(JL)*ZDIV65(JL)
  ENDDO
ELSEIF (KTILE  ==  4) THEN
  DO JL=KIDIA,KFDIA
    ZDIV75(JL) = 1.0_JPRB/(PWETL5(JL)+PRAQ5(JL))
    PQSAM5(JL) = PQS5(JL)+(PQMLEV5(JL)-PQS5(JL))*PWETL5(JL)*ZDIV75(JL)
  ENDDO
ELSEIF (KTILE == 6) THEN ! I.E. HIGH VEGETATION, SNOW FREE
  DO JL=KIDIA,KFDIA
    ZDIV85(JL) = 1.0_JPRB/(PWETH5(JL)+PRAQ5(JL))
    PQSAM5(JL) = PQS5(JL)+(PQMLEV5(JL)-PQS5(JL))*PWETH5(JL)*ZDIV85(JL)
  ENDDO
ELSE ! I.E. HIGH VEGETATION WITH SNOW (7)
  DO JL=KIDIA,KFDIA
    ZDIV95(JL) = 1.0_JPRB/(PWETHS5(JL)+PRAQ5(JL))
    PQSAM5(JL) = PQS5(JL)+(PQMLEV5(JL)-PQS5(JL))*PWETHS5(JL)*ZDIV95(JL)
  ENDDO
ENDIF

!*         2.5   DRY STATIC ENERGY AT THE SURFACE

DO JL=KIDIA,KFDIA
  PCPTS5(JL) = PTSKM1M5(JL)*RCPD*(1.0_JPRB+RVTMP2*PQSAM5(JL))
ENDDO

!     ------------------------------------------------------------------

!                 ----- ADJOINT COMPUTATION ------

!     ------------------------------------------------------------------

! Initialization

ZLIQ(:,:) = 0.0_JPRB

!*         2.5   DRY STATIC ENERGY AT THE SURFACE

DO JL=KIDIA,KFDIA
  PTSKM1M(JL) = PTSKM1M(JL)+RCPD*(1.0_JPRB+RVTMP2*PQSAM5(JL))*PCPTS(JL)
  PQSAM(JL) = PQSAM(JL)+RCPD*RVTMP2*PTSKM1M5(JL)*PCPTS(JL)
  PCPTS(JL) = 0.0_JPRB
ENDDO

!*         2.4   APPARENT SURFACE HUMIDITY

IF (KTILE  ==  1.OR. KTILE  ==  2.OR. KTILE  ==  3.OR. KTILE  ==  5) THEN
  DO JL=KIDIA,KFDIA
    PQS(JL) = PQS(JL)+PQSAM(JL)
    PQSAM(JL) = 0.0_JPRB
  ENDDO
ELSEIF (KTILE  ==  8) THEN
  DO JL=KIDIA,KFDIA
    PQS(JL) = PQS(JL)+PQSAM(JL)
    PWETB(JL) = PWETB(JL)+(PQMLEV5(JL)-PQS5(JL))*ZDIV65(JL)*PQSAM(JL)
    PQMLEV(JL) = PQMLEV(JL)+PWETB5(JL)*ZDIV65(JL)*PQSAM(JL)
    PQS(JL) = PQS(JL)-PWETB5(JL)*ZDIV65(JL)*PQSAM(JL)
    PWETB(JL) = PWETB(JL)-(PQMLEV5(JL)-PQS5(JL))*PWETB5(JL) &
     & *ZDIV65(JL)*ZDIV65(JL)*PQSAM(JL)
    PRAQ(JL) = PRAQ(JL)-(PQMLEV5(JL)-PQS5(JL))*PWETB5(JL) &
     & *ZDIV65(JL)*ZDIV65(JL)*PQSAM(JL)
    PQSAM(JL) = 0.0_JPRB
  ENDDO
ELSEIF (KTILE  ==  4) THEN
  DO JL=KIDIA,KFDIA
    PQS(JL) = PQS(JL)+PQSAM(JL)
    PWETL(JL) = PWETL(JL)+(PQMLEV5(JL)-PQS5(JL))*ZDIV75(JL)*PQSAM(JL)
    PQMLEV(JL) = PQMLEV(JL)+PWETL5(JL)*ZDIV75(JL)*PQSAM(JL)
    PQS(JL) = PQS(JL)-PWETL5(JL)*ZDIV75(JL)*PQSAM(JL)
    PWETL(JL) = PWETL(JL)-(PQMLEV5(JL)-PQS5(JL))*PWETL5(JL) &
     & *ZDIV75(JL)*ZDIV75(JL)*PQSAM(JL)
    PRAQ(JL) = PRAQ(JL)-(PQMLEV5(JL)-PQS5(JL))*PWETL5(JL) &
     & *ZDIV75(JL)*ZDIV75(JL)*PQSAM(JL)
    PQSAM(JL) = 0.0_JPRB
  ENDDO
ELSEIF (KTILE == 6) THEN ! I.E. HIGH VEGETATION, SNOW FREE
  DO JL=KIDIA,KFDIA
    PQS(JL) = PQS(JL)+PQSAM(JL)
    PWETH(JL) = PWETH(JL)+(PQMLEV5(JL)-PQS5(JL))*ZDIV85(JL)*PQSAM(JL)
    PQMLEV(JL) = PQMLEV(JL)+PWETH5(JL)*ZDIV85(JL)*PQSAM(JL)
    PQS(JL) = PQS(JL)-PWETH5(JL)*ZDIV85(JL)*PQSAM(JL)
    PWETH(JL) = PWETH(JL)-(PQMLEV5(JL)-PQS5(JL))*PWETH5(JL) &
     & *ZDIV85(JL)*ZDIV85(JL)*PQSAM(JL)
    PRAQ(JL) = PRAQ(JL)-(PQMLEV5(JL)-PQS5(JL))*PWETH5(JL) &
     & *ZDIV85(JL)*ZDIV85(JL)*PQSAM(JL)
    PQSAM(JL) = 0.0_JPRB
  ENDDO
ELSE ! I.E. HIGH VEGETATION WITH SNOW (7)
  DO JL=KIDIA,KFDIA
    PQS(JL) = PQS(JL)+PQSAM(JL)
    PWETHS(JL) = PWETHS(JL)+(PQMLEV5(JL)-PQS5(JL))*ZDIV95(JL)*PQSAM(JL)
    PQMLEV(JL) = PQMLEV(JL)+PWETHS5(JL)*ZDIV95(JL)*PQSAM(JL)
    PQS(JL) = PQS(JL)-PWETHS5(JL)*ZDIV95(JL)*PQSAM(JL)
    PWETHS(JL) = PWETHS(JL)-(PQMLEV5(JL)-PQS5(JL))*PWETHS5(JL) &
     & *ZDIV95(JL)*ZDIV95(JL)*PQSAM(JL)
    PRAQ(JL) = PRAQ(JL)-(PQMLEV5(JL)-PQS5(JL))*PWETHS5(JL) &
     & *ZDIV95(JL)*ZDIV95(JL)*PQSAM(JL)
    PQSAM(JL) = 0.0_JPRB
  ENDDO
ENDIF

!*         2.3   DEFINITION OF THE STOMATAL RESISTANCE AND BARE SOIL RES
!*               DOES WORK FOR TYPE 4, 6 AND 8 WHEN ROUTINE IS CALLED FOR 
!*               TYPE 4

IF (KTILE == 4) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV5(JL) > PQS5(JL)) THEN
      PWETL (JL) = 0.0_JPRB
    ENDIF
  ENDDO
ELSEIF (KTILE == 6) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV5(JL) > PQS5(JL)) THEN
      PWETH (JL) = 0.0_JPRB
    ENDIF
  ENDDO
ELSEIF (KTILE == 7) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV5(JL) > PQS5(JL)) THEN
      PWETHS (JL) = 0.0_JPRB
    ENDIF
  ENDDO
ELSEIF (KTILE == 8) THEN
  DO JL=KIDIA,KFDIA
    IF (PQMLEV5(JL) > PQS5(JL)) THEN
      PWETB (JL) = 0.0_JPRB
    ENDIF
  ENDDO
ENDIF

IF (KTILE  ==  4) THEN
  DO JL=KIDIA,KFDIA
    ZWROOTL = 0.0_JPRB
    ZWROOTH = 0.0_JPRB
    ZF2LP1 = 0.0_JPRB
    ZF21LP1 = 0.0_JPRB
    ZF2LP2 = 0.0_JPRB
    ZF21LP2 = 0.0_JPRB
    ZF2L   = 0.0_JPRB
    ZF21L   = 0.0_JPRB

    ZF2HP1 = 0.0_JPRB
    ZF2HP2 = 0.0_JPRB
    ZF2H   = 0.0_JPRB
    ZF21HP1 = 0.0_JPRB
    ZF21HP2 = 0.0_JPRB
    ZF21H   = 0.0_JPRB

    ZF2BP1 = 0.0_JPRB
    ZF2BP2 = 0.0_JPRB
    ZF2B   = 0.0_JPRB
    ZF21BP1 = 0.0_JPRB
    ZF21BP2 = 0.0_JPRB
    ZF21B   = 0.0_JPRB

    ZSRFL  = 0.0_JPRB
    ZF1LP  = 0.0_JPRB
    ZF1L   = 0.0_JPRB
    ZF1H   = 0.0_JPRB
    Z5S    = 0.0_JPRB
    Z6S    = 0.0_JPRB
    ZQSAIR = 0.0_JPRB
    ZCOR   = 0.0_JPRB
    ZF3L   = 0.0_JPRB
    ZF3H   = 0.0_JPRB
    ZWETL  = 0.0_JPRB
    ZWETH  = 0.0_JPRB

!           atmospheric moisture deficit stress function : F3

    ZF2B = ZF2B-RVRSMINB2D(JL)*ZF2BI5(JL)*ZF2BI5(JL)*PWETB(JL)
    PWETB(JL) = 0.0_JPRB
    PWETH(JL) = PWETH(JL)+PWETHS(JL)
    PWETHS(JL) = 0.0_JPRB
    IF (ZLAIH(JL) /= 0.0_JPRB) THEN
      ZWETH = ZWETH-ZRSMINH(JL)*ZWETHI5(JL)*ZWETHI5(JL)*PWETH(JL)
      ZF3H = ZF3H+ZLAIH(JL)*ZF1H5(JL)*ZF2H5(JL)*ZWETH
      ZF2H = ZF2H+ZLAIH(JL)*ZF1H5(JL)*ZF3H5(JL)*ZWETH
      ZF1H = ZF1H+ZLAIH(JL)*ZF2H5(JL)*ZF3H5(JL)*ZWETH
    ENDIF
    PWETH (JL) = 0.0_JPRB
    IF (ZLAIL(JL) /= 0.0_JPRB) THEN
      ZWETL = ZWETL-ZRSMINL(JL)*ZWETLI5(JL)*ZWETLI5(JL)*PWETL(JL)
      ZF3L = ZF3L+ZLAIL(JL)*ZF1L5(JL)*ZF2L5(JL)*ZWETL
      ZF2L = ZF2L+ZLAIL(JL)*ZF1L5(JL)*ZF3L5(JL)*ZWETL
      ZF1L = ZF1L+ZLAIL(JL)*ZF2L5(JL)*ZF3L5(JL)*ZWETL
    ENDIF
    PWETL(JL) = 0.0_JPRB
    IF (ZF3H55(JL) < ZEPSF3 .OR. ZF3H55(JL) > 1.0_JPRB) THEN
      ZF3H  = 0.0_JPRB
    ENDIF
    IF (ZF3L55(JL) < ZEPSF3 .OR. ZF3L55(JL) > 1.0_JPRB) THEN
      ZF3L  = 0.0_JPRB
    ENDIF
    ZQSAIR = ZQSAIR-ZF3H5(JL)*ZHSTRH(JL)*ZF3H
    PQMLEV(JL) = PQMLEV(JL)+ZF3H5(JL)*ZHSTRH(JL)*ZF3H
    ZQSAIR = ZQSAIR-ZF3L5(JL)*ZHSTRL(JL)*ZF3L
    PQMLEV(JL) = PQMLEV(JL)+ZF3L5(JL)*ZHSTRL(JL)*ZF3L

    ZCOR = ZCOR+ZQSAIR55(JL)*ZQSAIR
    ZQSAIR = ZCOR45(JL)*ZQSAIR

    ZQSAIR = ZQSAIR+RETV*ZCOR45(JL)*ZCOR45(JL)*ZCOR
    Z6S = Z6S+R2ES*PAPHMS5(JL)*ZDIV55(JL)*ZDIV55(JL)*ZQSAIR
    PAPHMS(JL) = PAPHMS(JL)-R2ES*Z6S5(JL)*ZDIV55(JL)*ZDIV55(JL)*ZQSAIR
    Z5S = Z5S+Z6S5(JL)*Z6S
    IF (PTMLEV5(JL) > RTT) THEN
      PTMLEV(JL) = PTMLEV(JL)+R3LES*(RTT-R4LES)*ZDIV45(JL)*ZDIV45(JL)*Z5S
    ELSE
      PTMLEV(JL) = PTMLEV(JL)+R3IES*(RTT-R4IES)*ZDIV45(JL)*ZDIV45(JL)*Z5S
    ENDIF

!           radiation stress function (proposed by Alan Betts): ZF1

    ZF1L = ZF1L+ZF1H
    IF (ZF1LP5(JL) > 1.0_JPRB) THEN
      ZF1LP = ZF1LP-ZF1L5(JL)*ZF1L5(JL)*ZF1L
    ENDIF
    ZSRFL = ZSRFL+RRSF1A*ZF1LPI5(JL) &
     & *(1._JPRB-(1.+ZSRFL5(JL))*ZF1LPI5(JL))*ZF1LP
    PSRFD(JL) = PSRFD(JL)+ZSRFL/RRSF1B

!           soil moisture stress function : F2

    IF (LEVGEN) THEN
      JS=KSOTY(JL)
      IF (JS >= 1) THEN
        ZQWEVAP=1._JPRB/(RWCAPM3D(JL,1_JPIM)-RWPWPM3D(JL,1))
      ELSE
        ZQWEVAP = 0.0_JPRB
      ENDIF
    ELSE
      ZQWEVAP=RQWEVAP
    ENDIF

    ZF21B = ZF21B + 2._JPRB*(1._JPRB-ZF21B5(JL))*ZF2B

    IF (ZF21BP25(JL) > RCEPSW) THEN
      ZF21BP2 = ZF21BP2+ZF21B
    ENDIF
    IF (ZF21BP15(JL) < 1.0_JPRB) THEN
      ZF21BP1 = ZF21BP1+ZF21BP2
    ENDIF
    ZLIQ(JL,1) = ZLIQ(JL,1)+ZQWEVAP*ZF21BP1

    IF (PTSAM1M5(JL,1) > RTT ) THEN
!   if first soil layer temperature is freezing then shutdown transpiration
!      ZF2H  = 2._JPRB*ZF21H-2._JPRB*ZF21H5*ZF21H

      ZF21H = ZF21H + 2._JPRB*(1._JPRB-ZF21H5(JL))*ZF2H

      IF (ZF21HP25(JL) > RCEPSW) THEN
        ZF21HP2 = ZF21HP2+ZF21H
      ENDIF
      IF (ZF21HP15(JL) < 1.0_JPRB) THEN
        ZF21HP1 = ZF21HP1+ZF21HP2
      ENDIF
      ZWROOTH = ZWROOTH+ZQWEVAP*ZF21HP1

      ZF21L = ZF21L + 2._JPRB*(1._JPRB-ZF21L5(JL))*ZF2L

      IF (ZF21LP25(JL) > RCEPSW) THEN
        ZF21LP2 = ZF21LP2+ZF21L
      ENDIF
      IF (ZF21LP15(JL) < 1.0_JPRB) THEN
        ZF21LP1 = ZF21LP1+ZF21LP2
      ENDIF
!    ENDIF

      ZWROOTL = ZWROOTL+ZQWEVAP*ZF21LP1
! Endif moved here from above as it seems to cover both L and H cases 
    ENDIF


    ZLIQ(JL,1) = ZLIQ(JL,1)+ZROOT1H(JL)*ZWROOTH
    ZLIQ(JL,2) = ZLIQ(JL,2)+ZROOT2H(JL)*ZWROOTH
    ZLIQ(JL,3) = ZLIQ(JL,3)+ZROOT3H(JL)*ZWROOTH
    ZLIQ(JL,4) = ZLIQ(JL,4)+ZROOT4H(JL)*ZWROOTH

    ZLIQ(JL,1) = ZLIQ(JL,1)+ZROOT1L(JL)*ZWROOTL
    ZLIQ(JL,2) = ZLIQ(JL,2)+ZROOT2L(JL)*ZWROOTL
    ZLIQ(JL,3) = ZLIQ(JL,3)+ZROOT3L(JL)*ZWROOTL
    ZLIQ(JL,4) = ZLIQ(JL,4)+ZROOT4L(JL)*ZWROOTL
  ENDDO

!                Compute first liquid fraction of soil water to 
!                be used later in stress functions
!          CONTRIBUTION TO APPARENT ENERGY, TAKING INTO ACCOUNT
!          FREEZING/MELTING OF SOIL WATER.

  DO JK=KLEVS,1,-1
    DO JL=KIDIA,KFDIA
      ZF = 0.0_JPRB
      ZLIQP1 = 0.0_JPRB
      ZLIQP2 = 0.0_JPRB

      IF (LEVGEN) THEN
        JS=KSOTY(JL)
        IF (ZLIQP25(JL,JK) > RWPWPM3D(JL,JK)) THEN
          ZLIQP2 = ZLIQP2+ZLIQ(JL,JK)
        ENDIF
        ZLIQ(JL,JK) = 0.0_JPRB

        IF (ZLIQP15(JL,JK) < RWCAPM3D(JL,JK)) THEN
          ZLIQP1 = ZLIQP1+ZLIQP2
        ENDIF
      ELSE
        IF (ZLIQP25(JL,JK) > RWPWP) THEN
          ZLIQP2 = ZLIQP2+ZLIQ(JL,JK)
        ENDIF
        ZLIQ(JL,JK) = 0.0_JPRB

        IF (ZLIQP15(JL,JK) < RWCAP) THEN
          ZLIQP1 = ZLIQP1+ZLIQP2
        ENDIF
      ENDIF
      PWSAM1M(JL,JK) = PWSAM1M(JL,JK)+(1.0_JPRB-ZF5(JL,JK))*ZLIQP1
      ZF = ZF-PWSAM1M5(JL,JK)*ZLIQP1

      IF (PTSAM1M5(JL,JK) < RTF1.AND.PTSAM1M5(JL,JK) > RTF2) THEN
        PTSAM1M(JL,JK) = PTSAM1M(JL,JK)-0.5_JPRB*RTF4 &
         & *COS(RTF4*(PTSAM1M5(JL,JK)-RTF3))*ZF
      ENDIF
    ENDDO
  ENDDO

ENDIF

!*         2.2   SATURATION PARAMETERS,

DO JL=KIDIA,KFDIA
  Z3S = 0.0_JPRB
  Z4S = 0.0_JPRB
  ZCOR = 0.0_JPRB

  IF (PTSKM1M5(JL) > RTT) THEN
    ZCOR = ZCOR+R5LES*PQS5(JL)*ZDIV35(JL)*ZDIV35(JL)*PDQS(JL)
    PQS(JL) = PQS(JL)+R5LES*ZCOR5(JL)*ZDIV35(JL)*ZDIV35(JL)*PDQS(JL)
    PTSKM1M(JL) = PTSKM1M(JL)-2.0_JPRB*R5LES*PQS5(JL)*ZCOR5(JL) &
     & *ZDIV35(JL)*ZDIV35(JL)*ZDIV35(JL)*PDQS(JL)
  ELSE
    ZCOR = ZCOR+R5IES*PQS5(JL)*ZDIV35(JL)*ZDIV35(JL)*PDQS(JL)
    PQS(JL) = PQS(JL)+R5IES*ZCOR5(JL)*ZDIV35(JL)*ZDIV35(JL)*PDQS(JL)
    PTSKM1M(JL) = PTSKM1M(JL)-2.0_JPRB*R5IES*PQS5(JL)*ZCOR5(JL) &
     & *ZDIV35(JL)*ZDIV35(JL)*ZDIV35(JL)*PDQS(JL)
  ENDIF
  PDQS(JL) = 0.0_JPRB

  ZCOR = ZCOR+ZPQS5(JL)*PQS(JL)
  PQS(JL) = ZCOR5(JL)*PQS(JL)
 
  PQS(JL) = PQS(JL)+RETV*ZCOR5(JL)*ZCOR5(JL)*ZCOR
  Z4S = Z4S+R2ES*PAPHMS5(JL)*ZDIV25(JL)*ZDIV25(JL)*PQS(JL)
  PAPHMS(JL) = PAPHMS(JL)-R2ES*Z4S5(JL)*ZDIV25(JL)*ZDIV25(JL)*PQS(JL)
  PQS(JL) = 0.0_JPRB

  Z3S = Z3S+Z4S5(JL)*Z4S
  IF (PTSKM1M5(JL) > RTT) THEN
    PTSKM1M (JL) = PTSKM1M (JL)+R3LES*(RTT-R4LES)*ZDIV15(JL) &
     & *ZDIV15(JL)*Z3S
  ELSE
    PTSKM1M (JL) = PTSKM1M (JL)+R3IES*(RTT-R4IES)*ZDIV15(JL) &
     & *ZDIV15(JL)*Z3S
  ENDIF

ENDDO

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('VSURFSAD_MOD:VSURFSAD',1,ZHOOK_HANDLE)
END SUBROUTINE VSURFSAD
END MODULE VSURFSAD_MOD
